07« 


VARIANCE  REDUCTION  IN  THE  SIMULATION  OP 
STOCHASTIC  ACTIVITY  NETWORKS 

George  S.  Fishman 

Technical  Report  Mo.  UNC/ORSA/TR-82/7 
January  1983 


Curriculum  in  Operations  Research 
and  Systems  Analysis 

University  of  North  Carolina  at  Chapel  Hill 


DTIC 

ELECTE1 
FEB  9  1983 1 


A 


7 

This  research  was  supported  by  the  Office  of  Naval  Research  under  contract 
N00014-?^6-C-0302.  Reproduction  in  whole  or  part  is  permitted  for  any  purpose 
of  the  United  States  government. 


This  document  has  been  approved 
for  public  release  and  sale;  its 
distribution  is  unlimited. 


Abstract 


This  paper  describes  a  Monte  Carlo  method  based  on  the  theory  of 
tjuaslrandom  points  for  estimating  the  distribution  functions  and  means  of 
network  completion  time  and  shortest  path  time  in  a  stochastic  activity 
network.  In  particular,  the  method  leads  to  estimators  whose  variances 
converge  faster  than  1/K  ,  where  K  denotes  the  number  of  replications 
collected  in  the  experiment.  The  paper  demonstrates  how  accuracy  diminishes 
for  a  given  K  with  Increasing  dimensionality  of  the  network  and  shows  how  a 
procedure  that  uses^a  cutset  of  the  network  together  with  convolution  can 
reduce  dimensionality  and  Increase  accuracy.  Two  examples  illustrate  the 
benefits  of  using  quasirandom  points  together  with  a  cutset  and  then 
convolution. 
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Introduction 


Activity  network  analysis  is  a  fairly  standard  tool  in  Operations 
Research.  When  arc  completion  times  in  such  a  network  are  deterministic, 
well-established  algorithms  exist  for  finding  its  network  completion  time  and 
shortest  path  time  in  a  computationally  efficient  manner.  For  example,  see 
Elmaghraby  (1977)  and  Wagner  (1975).  However,  when  arc  completion  times  are 
stochastic,  analysis  becomes  considerably  more  difficult,  even  for  relatively 
small  networks.  Here  the  objective  often  is  to  compute  the  distribution 
functions  (d.fs.)  and  means  of  network  completion  time  and  shortest  path 
time,  given  the  d.fs.  of  the  statistically  independent  individual  arc 
completion  times. 

Because  severe  difficulties  exist  in  deriving  analytical  solutions  to 
these  problems,  many  analysts  have  turned  to  Monte  Carlo  methods  to  derive 
approximate  solutions.  Van  Slyke  (1963)  proposes  the  use  of  importance 
sampling  to  gain  efficiency  in  estimating  the  characteristics  of  networks. 
Martin  (1965)  describes  how  one  can  increase  statistical  efficiency  when  the 
arc  passage  time  distributions  are  polynomials.  Burt  and  Garman  (1971) 
describe  how  conditional  sampling  can  reduce  dimensionality  and  thus  sampling 
variation.  Carrying  this  idea  one  step  further  Slgal,  Pritsker  and  Solberg 
(1979,  1980)  show  how  one  can  use  a  cutset  of  a  network  to  effect  similar 
reductions.  The  use  of  antithetic  variates  has  also  been  proposed  in  a  number 
of  studies,  Including  Cheng  (1980),  Sullivan,  Hayya  and  Schaul  (1982)  and 
Grant  (1982). 

Although  the  benefits  of  these  proposals  are  well  established, 
they  all  lead  to  estimators  with  variances  0(1 /K)  for  K  Independent 
replications.  The  present  paper  shows  how  to  increase  the  rate  of  convergence 
of  this  variance  using  quasirandom  points  and  illustrates  the  success  of  the 


approach  for  two  networks  of  Jf“8  and  N*18  arcs.  The  theory  of  quasirandom 
points  has  Its  origin  In  the  numerical  evaluation  of  multivariable  Integrals. 
The  theory  concerns  the  Identification  of  vector  sequences  that  when  used  In 
numerical  Integration  lead  to  approximations  whose  absolute  errors  converge  to 
zero  faster  than  If  Independent  random  vectors  were  used.  Since  one  can 
represent  the  desired  d.fs.  and  means  as  multivariable  Integrals,  the  appeal 
of  quasirandom  points  Is  apparent.  Moreover,  combining  the  use  of 
quasirandom  points  with  other  approaches  such  as  cutsets  and  convolution 
leads  to  even  greater  advantage,  as  we  demonstrate. 

Section  1  presents  the  stochastic  network  problem  In  detail.  Section  2 
describes  how  one  can  view  the  problem  as  one  of  multivariable  numerical 
Integration.  Section  3  describes  how  one  would  tackle  the  problem  by  crude 
Monte  Carlo  methods  and  also  describes  the  benefits  of  conditional  sampling. 
Section  4  shows  the  benefits  of  using  cutsets.  Section  5  then  Indicates  how 
cutsets  together  with  convolution  can  significantly  reduce  the  dimensionality 
of  the  problem. 

Section  6  presents  a  short  discussion  of  the  principal  concepts  of 
multivariable  numerical  integration  using  quasirandom  points  and  shows  the 
extent  to  trtilch  known  results  apply  to  the  problem  at  hand.  Section  7  lays 
out  an  experiment  whose  design  Is  used  In  Section  8  with  two  examples  to  show 
the  benefit  of  using  quasirandom  points,  especially  when  combined  with 
convolution  and  the  cutset  approach. 


1.  Definitions 


Consider  an  acyclic  directed  network  with  a  single  source,  single  sink,  N 
arcs  and  L  paths.  Let  Xi,...,Xn  >  the  passage  tines  for  arcs  1,...,N  , 
be  Independent  random  variables  where  X^  has  distribution  function  (d.f.)  Fi 
on  [0,«)  and  inverse  distribution  function  G^(u)  -  min(x:  Ff(x)  >  u  , 

0  <  u  <  1].  The  completion  time  of  path  m  is 

1-m  ”  aimXl  *  liel  *1  (1) 

m 

where  aim»l  if  arc  1  is  on  path  m  ,  aim“0  otherwise  and  Im  denotes 
the  set  of  arcs  on  path  m  . 

The  principal  purpose  of  this  study  is  to  characterize  the  network 
completion  time 

T*  -  max(Tlt 


and  the  shortest  path  time 

T*  -  min(Ti,...,TL)  .  (3) 

For  T*  ,  characterizations  include; 

a.  pr(T*  <  t)  0  <  t  <  • 

b.  E  T* 


c.  pr(m  is  the  longest  path)  m  -  1,...,L 

d.  pr(T*  <  t)  0  <  t  <  • 

e.  E  T* 

f.  pr(m  Is  the  shortest  path)  m  -  1 . L  . 

The  present  study  concerns  the  estimation  of  a,  b,  d  and  e  .  For  expository 
convenience  the  main  body  of  the  paper  concerns  the  estimation  of  pr(T*  <  t). 
The  Appendix  contains  details  for  the  estimation  of  the  remaining  quantities. 


(4) 


Let  | In|  denote  the  cardinality  of  Im  .  If 
N  “  Inr-1  I  ^m! 

then  no  paths  share  a  common  arc  and  T*  and  T*  have  d.fs. 

F*(t)  -  FIjB(t)  (5) 

and 

L 

F*(t)  ■  1  -  |1  -  FIn<t>l  0  <  t  <  •  , 

Fj  being  the  convolution  of  the  |lm|  arcs  on  path  m  .  As  an  example, 
m 

suppose  that  arc  times  are  exponential  with 

*i<0  -  1  -  e-Xtt  Xt  >  0  1  e  Ta  •  (6) 

Then  one  has  for  distinct  A* 

ft  (t)  -  l  -  l  — - — - 

m  lei*  "jerJ^Xi/Xj) 

j*i  (7) 

Usually  the  summation  in  (4)  exceeds  N  In  which  case  the  derivations  of 
F*  and  F*  are  more  complicated  than  in  (5).  In  practice,  even  when  (4) 
holds,  the  convolutions  for  Fj^,...,Fj^  may  prove  too  complex  to  derive 
analytically. 

2.  A  Solution  via  Integration 

An  alternative  approach  to  characterization  Is  through  Integration;  In 

particular,  numerical  Integration.  Observe  that  one  can  write 

F*(t)  -  E  I[0,t)(T*) 

11  H 

"  ^  *  * *0  ^ ^**1 » • • • »*n^  ^jbj^7i(xi) 


(8a) 


M«,b)  (*)  -  1  if  a  <  x  <  b 

”  0  otherwise, 

N  N 

h(t;x! . *N>  "  I[0,t)  (“axCE^ataXi . J1Hl  «iL*i))  (9a) 

and 

g(t;«i, •  •  •  ,ujj)  -  h(t;G!(ui) . GN(uN))  .  (9b) 


Now  (8b)  Involves  multiple  integration  over  the  N-dimensional  unit  hypercube 
[0,l)x...x[0,l)  ,  hereafter  denoted  by  .  it  is  the  expression  (9b)  on 
which  our  analysis  is  based. 

Consider  the  approximation 

"  K  ^ j-1  8(t’uij»*,,*u||j)  (10) 

where  uj  -  (nij»  j  —  1,2,...}  i-l,...,N  are  infinite  sequences  of  points 
in  chosen  in  one  of  several  ways.  One  way  is  through  pure  random  sampling 

which  we  call  the  crude  Monte  Carlo  method . 

3.  Crude  Monte  Carlo  Sampling 
Let 

T*  -  max  (Tj1.-**,T;jL) 

where 

^Jm  "  Ii-1  aim  ®i(uij)  ®  “  1,...,L 

and  uij  i  -  1,...,N  j-1 . K  are  i.i.d.  from  U[0,1) 

8j(t)  -  g(t;ulj,...,uNj)  -  I[0,t)(Tj)  • 

Then  g^t)  in  (10)  has  expectation  F*(t)  and 
var  ,  (t)  -  P*(t)1  . 


(ID 

(12) 

.  Also  define 
(13) 
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To  Improve  on  the  result  In  (14) ,  Burt  and  Garman  (1971)  Introduce  the 

use  of  conditional  sampling.  Suppose  that  there  exist  arcs  li . ll,  unique 

to  paths  1,...,L  respectively.  Assume  that  these  arcs  are  distinct  and 


define 


sjm  *  alm  Ol(«lj) 

±*U 


8j(t)  “  nm-l  Fl“(t  _  .  (16) 

Now  g^(t)  in  (10),  based  on  (16),  has  the  correct  expectation  P*(t)  ,  as 

can  be  seen  by  Integrating  (16)  with  respect  to  the  joint  p.d.f.  of 

Sjp...,SjL  .  But  gR(t)  has  smaller  variance  than  (14),  as  a  result  of 

sampling  from  N-L  instead  of  N  arcs.  Also  note  that 

L 

s j(0  <  Plm(t)  (17) 

provides  an  upper  bound  which  also  applies  to  5g(t). 

In  principle,  conditional  sampling  can  be  made  more  effective  by  working 


with  subsets  of  arcs  unique  to  each  arc.  Let 
Jm  "  Tj) 

so  that  Jim  Is  the  set  of  arcs  unique  to  path  m 
m  ”  1,...,L  and  define 


m  *  1 , • . . ,L  (1? 
.  Assume  |Jm|  >  1 


sJm  “  <19> 

and 

L 

8j(t)  -  \ml  P^t-Sja)  (20) 

where  Ft  denotes  the  d.f.  of  the  sum  of  the  passage  times  on  the  arcs  unique 
m 

to  path  m  .  Since  sampling  in  N  +  L  -  l  ^ | -k|  dimensions  is  all  that  Is 
necessary,  the  use  of  (20)  in  (10)  preserves  unbiasedness  and  reduces  variance 


even  further. 


In  essence,  the  effect  of  conditional  sampling  Is  to  replace  the 
coefficient  of  1/K  in  (14)  by  one  of  smaller  magnitude  but  In  no  way  to 
change  the  convergence  rate  with  regard  to  K  as  K  ♦  «•  .  Also,  note  the 
upper  bound 

<  Vi  Vl)  <21) 

which  Is  tighter  than  the  bound  In  (17). 

4.  Cutsets 

In  practice  It  Is  conceivable  that  at  least  one  path  does  not  have  a 
unique  arc.  That  Is,  |Jm|  -  0  for  at  least  one  path  m  .  Figure  1 
Illustrates  such  a  case.  Although  conditional  sampling  as  described  by  Burt 


Insert  Fig.  1  about  here. 


and  Garraan  does  not  apply  here,  another  more  general  approach  proposed  by 
Slgal,  Pritsker  and  Solberg  (1979,  1980)  does.  Let  H  denote  a  cutset  of  the 
network.  A  cutset  Is  a  set  of  arcs  that  connects  a  set  of  nodes  W  containing 
the  source  with  a  set  of  nodes  (V  containing  the  sink.  Also,  assume  that 
each  path  has  only  one  arc  In  H  .  If  each  arc  In  H  points  from  W  to 
W  ,  H  is  called  a  uniformly  directed  cutset  (UDC) . 


Define 


r" 

sjm  “  li-1  ®lm  Gi(uij) 

1 

Yij  "  ®“P  (®Jm  Sjm) 

m-1  |  #  e  e  |L 


1  e  H  . 


Then  define 


81(t>  "  *i(t  -  Yi  )  • 


g_(t)  in  (10)  basea  a  4)  gives  an  unbiased  estimate  of  F*(t) 


with  smaller  variance  than  crude  Monte  Carlo  sampling  since  one  is  sampling  in 
only  N  -|W|  dimensions.  Since  a  network  may  have  more  than  one  UDC,  using 
the  one  with  maximal  cardinality  gives  the  greatest  reduction  in 
dimensionality.  However,  the  determination  of  this  cutset  for  an  arbitrarily 
large  network  is  not  trivial,  being  an  NP  complete  problem.  Nevertheless, 
among  alternative  known  cutsets  for  a  given  network,  the  one  with  largest 
cardinality  serves  our  purpose  best.  As  an  example  of  the  use  of  cutsets  note 
that  (e2,  e3,  e4,  ej)  and  (e3,  e5,  eg,  ej)  are  UDCs  with  maximal 
cardinality  4  for  the  network  of  Fig.  1,  whereas  the  cutset  (ej,  e2,  e3) 
with  cardinality  3  is  not. 

As  a  second  example,  consider  the  network  in  Fig.  2,  taken  from  Battersby 
(1970),  which  describes  the  steps  encountered  in  the  partial  overhaul  of  a 
unit  in  an  oil  refinery.  Table  1  shows  the  incidence  matrix  a  ■  | |  a£m  || 

Insert  Fig.  2  about  here. 

for  arcs  (rows)  and  paths  (columns).  Observe  that  arcs  (e2,  e3,  e5,  eg,  eg, 
e10»  e12)  and  arcs  (e2,  e3,  es,  eg,  ej.o>  ei2*  eis)  form  two  UDCs  with  maximal 
cardinality  7  .  However,  the  cutset  (e2,  e3,  e5>  eg  +  ei5,  eg,  e^o*  e12)  » 
formed  by  combining  arcs  eg  and  ei5  (see  Fig.  2)  has  cardinality  8  and 
can  be  used  to  advantage  to  reduce  dimensionality. 


Insert  Table  1  about  here. 


5.  Further  Dimensional  Reduction 

One  can  achieve  at  least  one  additional  reduction  in  the  dimensionality 
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J  ■  nr  t  -  yk 

ml . mr  j-1  ®J  Am-1 

*  •  •  •  »®r 


of  sampling  by  exploiting  available  knowledge  about  a  network.  Let 

Im  -  Io  "  Im  H  H  (25) 

where  H  denotes  a  UDC.  Then  define  the  subsets 

icn  ng  (26) 

1  <  mi  <  m2  <  ...  <  mr  <  L ;  r*l , . . . ,L  . 

Here  Jm  m  denotes  the  subset  of  arcs,  not  In  H  ,  that  are  uniquely  on 

1 . r 

paths  .  Suppose  there  are  q  such  nonempty  sets  Ri . Rq. 

Let  Fg  denote  the  d.f.  of  the  sum  of  the  |R^|  arc  passage  times 
corresponding  to  the  arcs  in  R^  and  let  Gg  denote  the  corresponding 
Inverse  d.f.  Then  define 

Sjm  *  (27) 

and 


Ytj  -  sup  (Sjm)  i-1 . |tf|  (28) 

m=*l . L 

and 

8j(t)  "  Vf/  Fi<c  -  TiJ>  •  (29) 

Now,  gR(t)  *n  (10)  based  on  (29)  again  is  an  unbiased  estimator  of 

V  3 

F*(t).  However,  since  N  -  |w|  +  q  -  |Rj|  dimensions  arise  for  sampling 

per  replication,  one  anticipates  that  var  gg(t)  has  smaller  magnitude  than 
previously  described  methods  produce.  Note  that  this  reduced  dimensionality 
comes  from  the  use  of  UDCs  together  with  convolution. 

6.  Multivariable  Numerical  Integration 

As  (8b)  shows,  one  can  view  F*(t)  as  the  result  of  a  multivariable 
integration  over  a  restricted  region  in  .  We  now  describe  how  taking  this 


view  of  the  problem  leads  to  useful  results  when  quasirandom  points  are  used. 
Our  account  follows  Kulpers  and  Hiederreiter  (1974)  and  Nlederrelter  (1978). 
Also  see  Schmidt  (1973). 

Suppose  that  our  objective  Is  to  evaluate  the  Integral 
1  1 


B 


/.../ 
0  0 


f(x 


I’*** 


■V 


dx  . . .dx 
1  N 


(30) 


by 

®K  m  if  I  (31) 

and  to  derive  bounds  on  the  error 

4C  -  |B  -  Br|  •  (32) 

Three  definitions  facilitate  our  description.  Consider  the  sequences 
tuj  ■  (u^j, • . . >u^j)i  j  m  1» • *  «  »K)  • 

Definition  1.  A(R;K)  -  number  of  points  «y.,  •••,»&  that  R  C  . 

Definition  2.  u^ . u^  are  uniformly  distributed  In  Rn  if 

11m  »  i 

K*»  KX( R) 

for  all  R-  {(xt,...,xN):  ai  <  xj  <  Si  1-1 . N)  R  C  and  where 

A(R)  Is  the  measure  of  volume  of  R  In  Vjj  . 

Definition  3.  The  extreme  discrepancy  associated  with  a  sequence  . . . 

Is 

dK  -  DK(jii . Wk)  "  jup  _  x(R*)| 

r*ci/n 

where  R*  -  {(Xj . xN):  0  <  X1  <  St  1  -  1, ...  ,19}  .  Hereafter  we  assume 

that  {iy  j  -  1,...,K)  is  a  uniformly  distributed  sequence  on  Rtf  . 

These  definitions  together  with  several  fundamental  theorems  from  the 


theory  of  quasirandom  points  provide  us  with  a  way  of  characterizing  the 
error  (32).  In  particular, 
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Theorem  1  (Hlawka  1961).  If  f  is  of  bounded  variation  in  the  sense  of 
Hardy  and  Krause,  then 

*K  <  V(f,N)DK  (33) 

where  V(f,N)  is  a  function  of  the  bounded  variation  in  f  in  N  and  lower 
dimensions. 

Theorem  2  (van  Aardenne-Ehrenfest  1945).  For  any  infinite  sequence 
with  N  >  1 


lim  KDr  m  «a  . 

K-h# 

Theorem  3  (Roth  1954)  .  Por  any  sequence  of  K  points  in  with  N  >  2 

C  (log  K)(N-l)/2 

dK  >  a— -g - 

where  Cfl  is  a  function  of  N  only. 

Theorem  4  (Roth  1954)  .  For  any  infinite  sequence  in  Vjj  with  N  >  1 
Cjj  (log  K)»/2 

”>=> - g - 

where  Ctf  is  a  function  of  N  only. 

Theorem  1  shows  that  one  can  bound  the  error  (32)  by  a  quantity 
proportional  to  the  extreme  discrepancy  Dr  ,  a  most  convenient  result. 
Theorem  2  implies  that  one  cannot  expect  the  extreme  discrepancy  to  converge 
as  fast  as  1/K  .  Theorems  3  and  4  provide  lower  bounds  on  how  fast 
convergence  can  occur. 

We  now  describe  a  uniformly  distributed  sequence  for  which  upper  bounds 
are  known.  If  R  >  2  is  an  integer,  then  every  non-negative  integer  n  has 
an  expansion  of  the  form 


ai  e  (0,1,..., R-l)  (34) 

0  <  i  <  m  and  m  -  LlogRnJ  . 


n  "  IT-o  aiRi 
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Corresponding  to  this  expansion  one  has  the  radical  Inverse  function 

*R(n)  -  I“,0  aiR-1'1  .  (35) 

Ue  now  define  several  sequences  based  on  (34): 

van  der  Cor put  Sequence:  t+2 ( j) ;  j  ■  0,1,...,K-1} 

Haaaersley  Sequence:  {^(j),...,  to^U),  i  j  •  0,1, . . .,K-1 } 
where  Rl,...,RR-l  are  pairwise  relatively  prime. 

Hal  ton  Sequence:  (^(J),...,  j);  j-0,1 . K-l  } 

where  Ri,...,Rfl  are  pairwise  relatively  prlae. 

The  van  der  Corput  sequence  Is  unlforaly  distributed  on  and  was  the  first 

sequence  of  this  kind  proposed  for  univariate  nuaerlcal  Integration.  The 
Hammersley  and  Halton  sequences  are  unlforaly  distributed  on  and  as 

Theoreas  5  and  6  show  they  offer  useful  upper  bounds. 

Theorem  5  (Halton  1960).  For  the  Haaaersley  sequence 


.  (log  K)N-1  Jf-1  ,3Ri-2  v 
Dk  <  — —  n‘1.1  (I5rSl) 

Theorea  6  (Halton  1960).  For  the  Halton  sequence 

(log  K)«  JH1  *1^2 
^  <  K  “l-l  ' log  Ri} 


N  >  max  Ri 


N  >  ajx  Ri  . 


Several  Issues  of  significance  need  to  be  mentioned  now.  Firstly,  since 
the  van  der  Corput,  Haaaersley  and  Halton  sequences  are  deterministic  the 
upper  bounds  can  be  Interpreted  as  worst  case  bounds  with  certainty. 

Secondly,  although  V(f,N)  In  (33)  Is  In  theory  computable,  In  the  present 
case  f  Is  unknown.  Thirdly,  the  bounds  presented  here  hold  for  Integration 
the  entire  H-dlaenslonal  hypercube;  however,  as  (8a),  (9a)  and  (9b)  show,  our 
problem  calls  for  Integration  over  a  restricted  region  In  .  Now  It  Is 

known  (Nlederreiter  1978,  p.  982)  that  for  integration  over  an  arbitrary 
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region  In  Vn 

*  <  C''(Dk)1/h  (36) 

where  CjJ'  Is  a  function  of  N  .  This  is  a  considerably  weaker  bound  but  since 
it  applies  to  such  a  wider  set  of  regions  in  than  (33)  does  there  is 

reason  to  believe  that  the  use  of  the  Hamersley  and  Halton  sequences  can  lead 
to  accelerated  convergence. 

The  fourth  issue  concerns  boundedness.  In  the  case  of  estimating  F*(t) 
and  F*(t)  ,  this  assumption  is  met.  However,  it  is  not  met  for  ET*  and 
ET*  .  This  fact  together  with  the  integration  over  a  restricted  region  in  Vtf 
makes  clear  that  when  we  move  from  theory  to  practice,  some  skepticism  exists 
as  to  how  efficient  our  approach  will  be.  As  the  examples  in  Section  8  show, 
there  is  little  basis  for  this  skepticism. 

The  fifth  issue  relates  to  dimensionality.  As  Theorems  5  and  6  make 
clear,  it  is  in  one's  best  Interest  to  reduce  dimensionality  as  much  as 
possible  before  using  quasirandom  points.  The  cutset  approach  together 
with  convolution  are  the  methods  we  use  to  effect  this  reduction. 

The  sixth  issue  of  significance  concerns  the  choice  between  the 
Hammersley  (finite)  and  Halton  (infinite)  sequences  in  practice.  If  one  uses 
the  Hammersley  sequence  for  fixed  K  and  then  decides  that  the  accuracy  is 
not  acceptable,  there  is  no  recourse  to  continued  exploitation  of  its  special 

structure  since  0,1/K . (K-l)/K  are  the  values  in  the  Nth  dimension.  By 

contrast,  the  Halton  sequence  enables  us  to  continue  exploitation  by  merely 

generating  additional  points  {^(j) . ^(j)}  j-k,  k+1 .  The  modest 

degradation  in  the  upper  bound  that  arises  when  using  the  Halton  rather  than 


the  Hammersley  sequence  seems  worth  the  price. 

7.  Experimental  Design 

This  section  describes  the  layout  for  an  experiment  designed  to  determine 
the  extent  to  which  quasirandom  points  lead  to  accelerated  accuracy 
when  estimating  F*(t)  ,  ET*  ,  F*(t)  and  ET*  for  the  networks  in  Figs.  1 
and  2  .  In  particular,  we  introduce  a  degree  of  randomness  into  the 
experiment  in  order  to  compute  estimates  of  the  variances  of  our  point 
estimates.  We  then  study  the  behavior  of  these  sample  variances  as  K 
increases . 

f  nsider  an  experiment  conalsting  of  Q  statistically  Independent  blocks 
or  macroreplications  each  of  K  microreplications .  Let  {Ula:  i-l,...,N; 
m-1, . . . ,Q}  denote  a  sequence  of  i.i.d.  random  variables  from  u[0,l)  and 
define 

sim  -  (s:  ♦RjU)  -  Ui«)  i-l,...,N  hf1,...,Q  .  (37) 

Then  on  macroreplication  m  we  use  the  quasirandom  point 

♦^(Sim  +  j-1 ),..., <t»RN(SMn  +  j-1)  on  microreplication  J  j«l,...,K  . 

Let  8jm  denote  the  unbiased  estimate  of  a  particular  quantity  0 
computed  on  microreplication  j  on  macroreplication  m  .  Then  01,..*,^, 
where 

Ao  ■  K  m*l,...,Q  ,  (38) 

are  i.i.d.  random  variables  with  sample  variance 

\  -  <Fi  Cl  (5.  •  V2  <”> 

where 

\  -  $  ^  ■ 

It  is  the  behavior  of  s£  versus  K  that  interests  us. 


To  provide  an  Instructive  normalization  we  also  run  M  statistically 


Independent  macrorepllcatlons  each  of  K-l  mlcrorepllcatlon.  Let  8^ . ^ 


denote  the  resulting  estimates  of  6  each  with  sample  variance 

wq  “  m=t  ^2-i  ( V  V2 

where 

1  cM 


(40) 


8  .if  a 

m  M  4*-l  ^  * 

Then  the  quantity  w^  /Ks£  should  Increase  as  K  Increases  If  accelerated 


convergence  Is  occurring. 

We  set  M-10^  .  For  the  quasirandom  case  we  perform  runs  for  K-2^  , 
j-1,2,...  and  set  Q-216/K  for  K<29  and  Q-2?  for  K>29  .  Following  this 


procedure  enables  us  to  treat  s£  and  w£  as  highly  reliable  point 


estimators  of  the  quantities  var  8*  and  var  6^  respectively. 

8.  Examples 

We  refer  to  the  network  In  Fig.  1  as  the  SPS  network  and  to  the  one  In 
Fig.  2  as  the  Battersby  network.  Table  2  lists  the  relevant  cutsets  and 
noncutsets.  In  both  examples  the  arc  passage  times  are  assumed  Independent 


Insert  Table  2  about  here. 


and  exponentially  distributed.  For  the  SPS  network  the  means  are  Xi-10  , 
*2"15  ,  \3«18  ,  X4«6  ,  Xs«7  ,  X6-3  ,  X7-IO  and  X8"2  .  For  arc  1  on 
mlcrorepllcatlon  j  on  macroreplication  m  the  arc  passage  time  Is 

*1  "  GiCfc^Sta  +  J-l))  -  -Xiln(l  -  ^(Si,  +  j-l»  .  (41) 

For  the  runs  using  quasirandom  points,  the  cutset  W-(e2,  e3,  e4,  e7)  was 
used  reducing  dimensionality  to  N-|#|  ■  8-4-4  .  An  algorithm  due  to  Halton 


and  Smith  (1964)  was  used  to  coapute  the  Halton  sequences.  For  the  M 
independent  aicroreplications  all  N  arc  times  froa  (41)  were  uaed. 

Table  3  shows  the  behavior  of  w2/Ks^  for  F*(t)  estimated  at 
t-15,  30,  40,  50,  70  for  ET*  ,  for  F*(t)  estimated  at  t-5,  10,  15,  20,  25 
and  for  ET*  .  As  a  point  for  coaparison,  Slgal,  Prltsker  and  Sol berg  (1979) 

Insert  Table  3  about  here. 

report  variance  reductions  ranging  from  4.18  to  11.00  for  the  estimation 
of  F*(t)  at  t*15,  30,  40,  50  and  70  using  the  cutset  ff  and  independent 
aicroreplications . 

The  second  example  studies  variance  reduction  for  the  Battersby  network 
using  the  means  in  Table  1.  TWo  separate  sets  of  experiments  were  run.  The 
first  set  uses  partition  1  in  Table  2  with  the  cutset  H~(e2»  e3»  e5»  e6» 
e9,  eio>  ei2)  with  quasirandom  points,  yielding  a  dimensionality  of 
N-|tf|~ll  .  The  second  set  uses  partition  2  with  the  cutset  Hm(e2>  e3»  e5» 
e6  +  eis,  eg,  eio,  ei2)  yielding  a  dimensional ity  of  6  after  appropriate 
convolution  within  the  noncutset.  The  quantity  F*(t)  was  estimated  for 
t-90,  110,  120,  130,  140,  150,  160,  180,  190  and  F*(t)  for  t-20,  25,  30, 
35,  40,  50,  60,  80,  90  .  Table  4  shows  the  resulting  variance  reductions 
wg/Ks£  for  F*(t)  ,  ET*  ,  F*(t)  and  ET*  . 


Insert  Table  4  about  here. 
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Tables  3  and  4  adalt  several  instructive  observations: 

1.  As  K  Increases  variance  reduction  Increases;  however, 
variance  reduction  is  not  necessarily  aonotone  in  K. 

4.  For  given  K  ,  variance  reduction  tends  to  decline  with 
increasing  dimensionality. 

3.  For  given  K  variance  reduction  is  greater  for  F*(t)  than 
for  ET*  and  for  F*(t)  than  for  ET*  ,  showing  the  benefit 
one  obtains  froa  bounded  integrands. 

Although  these  observations  establish  the  benefits  of  using  quaslrandoa 
points  with  the  cutset  approach  and  convolution,  there  is  a  cost  issue  that 
also  needs  to  be  considered.  If  one  coapares  the  cost  of  producing  one  saaple 
point  for  F*(t)  ,  ET*  ,  F*(t)  and  ET*  using  quaslrandoa  points  to  the 
cost  of  producing  a  saaple  point  using  pure  randoa  saapllng,  the  ratio  never 
exceeds  2.5  in  our  experlaents.  The  aost  costly  step  turns  out  to  be 
saapllng  the  suamed  arc  tiaes  xj  +  Xj  t  in  the  noncutset  of  partition  2  in 
Table  2,  froa  the  distribution  function 

-X4t  -X.it 

r<t>  •  1  -  FQij  '  PxyxT  ' 

This  was  done  using  the  Newton-Raphson  aethod. 

Since  this  cost  ratio  is  Independent  of  K  ,  we  have  strong  evidence 
that  for  sufficiently  large  K  ,  the  quaslrandoa  point  aethod  offers  a  clear 
advantage  over  pure  random  saapllng  for  networks  of  arbitrary  size.  A  prograa 
called  NETWOK  is  currently  under  developaent  to  lapleaent  the  proposed  aethod 
in  a  computationally  efficient  manner. 
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Append  lx 

Eatlmatlon  of  F*(t).  We  derive  the  estimator  based  on  the  cutset  H  and 
(22).  Let 

Zij  -  inf  (S ja)  1  e  H 

m*l |  •  *  •  jL 

aim-i 

and  define 

hj(t)  -  1  -  nieH  [1  -  Pi(t  -  Zij)]  . 

Then 

*»K(t)  -  ^1^.1  *>j(t) 

gives  an  unbiased  estimator  of  F*(t)  based  on  sampling  from  N  -  |  ff| 
dimensions. 

Estimation  of  ET*  and  ET*.  Using  (11)  one  has  the  unbiased  estimator 
S*  1  rK  _* 

tk  -  TT  Xj-i  fj 
for  ET*  .  Let 

T*j  -  min(Tji,...,TjL)  • 

Then 

T*K  -  IT  I  j-l  T*  j 

is  an  unbiased  estimator  of  ET*  . 


Table  2 


Network 


Fig.  1 


Battersby  (Fig.  2) 
Partition  1 


Battersby  (Fig.  2) 
Partition  2 


ei  +  ej  =  convolution 


Network  Analysis 

Cutset  Noncutset 


c2,  e3,  e4,  e7 


•1*  e5*  eg*  e8 


e2»  e3,  es,  e6,  eg,  eio, 

e12 


ei,  e4,  e7,  es,  en 
ei3»  e14»  eis,  ei6, 
e17 •  eis 


e2»  ®3»  es,  e6  +  eis,  eg, 

«10.  ei2 


ei  +  e4,  e7  +  eifj, 
e8  +  en,  ei3,  ei4, 
ei7  +  ei8 


of  the  d.fs.  of  arcs  i  and  j  . 


Variance  Reduction  w2  /Ks£  for  the  SPS  Networkt 


K 

F* 

min 

M-105  , 

(t) 

max 

N-8  ,  L- 

ET* 

■6  ,  |h|-4 

F*(t) 

min  max 

ET* 

2 

42 

106 

7 

46 

196 

8 

22 

53 

96 

9 

44 

345 

9 

23 

74 

101 

12 

53 

773 

11 

2* 

87 

177 

17 

125 

957 

19 

25 

111 

258 

23 

149 

1,285 

24 

26 

146 

701 

27 

313 

2,722 

36 

27 

248 

1,066 

51 

534 

2,984 

45 

28 

365 

1,409 

83 

710 

3,183 

51 

29 

383 

2,426 

130 

1,527 

5,323 

123 

987 

4,612 

153 

2,197 

4,512 

138 

2H 

1,000 

5,221 

266 

3,438 

8,604 

260 

212 

1,518 

11,341 

338 

5,391 

13,213 

303 

213 

1,978 

16,907 

548 

17,299 

39,357 

484 

214 

3,289 

51,709 

1,258 

27,943 

63,234 

838 

^mln  and  max  for  F*(t)  denote,  respectively,  the  minimal  and  maximal 
variance  reduction  for  t*15,  30,  40,  50,  70  .  min  and  max  for  F*(t) 
denote,  respectively,  the  minimal  and  maximal  variance  reduction  for 
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maximal  variance  reduction 
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